home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / lisp / eulisp / mpfeel.lha / MPFeel / Plurals / Random / fp_random.m < prev    next >
Text File  |  1992-05-12  |  2KB  |  70 lines
  1. /******************************************************************************
  2.  
  3.    A portable plural random number generator based on three linear
  4.    congruential generators.  -- Numerical Recipes in C pp. 210
  5.  
  6.    Copyright 1991,1992 Lee Iverson
  7.  
  8.    Author:
  9.      Lee Iverson <leei@mcrcim.mcgill.edu>,
  10.      McGill Research Centre for Intelligent Machines (McRCIM)
  11.  
  12.    See "copyright.h" for complete copyright information.
  13.  
  14. ******************************************************************************/
  15.  
  16. static char RCSid[] = 
  17.   "$Id: fp_random.m,v 1.2 92/01/23 17:31:47 leei Exp $";
  18.  
  19. #include <values.h>
  20. #include <math.h>
  21.  
  22. #include "p_randomI.h"
  23.  
  24. plural float fp_frandom _IMPL_VOID()
  25. {
  26.     register plural float tmp;
  27.     register plural int j;
  28.  
  29.     if ( !rand_init__ ) init_p_random(1 /* Was 0 */);
  30.  
  31.     ix1__ = (IA1*ix1__ + IC1) % M1;
  32.     ix2__ = (IA2*ix2__ + IC2) % M2;
  33.     ix3__ = (IA3*ix3__ + IC3) % M3;
  34.     j = 1 + ((97*ix3__)/M3);
  35.     tmp = r__[j];
  36.     r__[j] = (ix1__+ix2__*RM2)*RM1;
  37.     return tmp;
  38. }
  39.  
  40. plural float
  41. fp_uniform_random _IMPL((lo, hi),
  42.             register plural float lo _AND
  43.             register plural float hi)
  44. {
  45.     return lo + fp_frandom()*(hi-lo);
  46. }
  47.  
  48. plural float
  49. fp_normal_random _IMPL((stdev), plural float stdev)
  50. {
  51.     static int iset = 0;
  52.     static plural float gset;
  53.     plural float fac, r, v1, v2;
  54.  
  55.     if ( !iset ) {
  56.     do {
  57.         v1 = fp_uniform_random((plural float)-1.0,(plural float)1.0);
  58.         v2 = fp_uniform_random((plural float)-1.0,(plural float)1.0);
  59.         r = v1*v1 + v2*v2;
  60.     } while ( r >= 1.0 );
  61.     fac = fp_sqrt(-2.0*fp_log(r)/r);
  62.     gset = v1*fac; ++iset;
  63.     return v2*fac*stdev;
  64.     } else {
  65.     --iset;
  66.     return gset*stdev;
  67.     }
  68. }
  69.  
  70.